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FIELD OF THE INVENTION 

The present invention relates to a method of forming, from data obtained by 
exploration of a zone of a heterogeneous medium, a model representative of the 
distribution in the zone of a physical quantity (at least partly) free of the presence of 
5 correlated noises that may be contained in the data. 

The method applies for example to the quantification of the acoustic impedance in 
an underground zone. 

BACKGROUND OF THE INVENTION 

The process consisting in seeking a model that adjusts to experimental 
10 measurements has developed in nearly all the fields of the sciences or technology. Such 
an approach is known under various names : least-squares method for parameters 
estimation, for inverse problem solution. For a good presentation of this approach 
within the context of geosciences, one may for example refer to : 

- Tarantola, A. : « Inverse Problem Theory : Method for Data Fitting and Model 
15 Parameter Estimation », Elsevier, Amsterdam, 1987. 

It can be noted that the term « lest squares » refers to to the square of the norm in 
the data space for quantifying the difference between the response of a model (which is 
the image of the model by a previously selected modelling operator) and the data, a cost 
function that has to be minimized to solve the problem. Using the square of the norm to 
20 define the cost function is just a practical convenience, and it is not fundamentally 
essential. Besides, many authors use, for various reasons, another definition of the cost 
function but this definition remains based on the use of the norm, or of a semi-norm, in 
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the data space. Finally, we have considerable latitude for selecting the norm (or the 
semi-norm) in the data space (we are in no way forced to use the Euclidean norm). In 
the case of noise-containing data, the solution can substantially depend on the choice 
made at this stage. For more developments concerning this problem, one can for 
5 example refer to the following publications : 

- Tarantola, A. : « Inverse Problem Theory : Method for Data Fitting and Model 
Parameter Estimation », Elsevier, Amsterdam, 1987 ; Renard and Lailly, 2001 ; Scales 
and Gersztenkorn, 1988 ; Al-Chalabu, 1992. 

The measuring results often contain errors. Modelling noises add further to these 
10 measuring errors when the experimentation is compared with modelling results : 
modellings are never perfect and therefore always correspond to a simplified view of 
reality. The noise will therefore be described hereafter as the superposition of : 

- non correlated components (white noise for example), 

correlated components, which means that the existence of a noise on a measuring 
15 sample translates into the existence of a noise of equal nature on certain neighbouring 
measuring points ; modelling noises typically belong to this category. 

When the data contain correlated noises, the quality of the model estimated by 
solving the inverse problem can be seriously affected thereby. As already mentioned, no 
modelling operator is perfect. It is therefore the work of the whole community of the 
20 people involved in the identification of parameters describing a model which is 
hindered by the existence of correlated noises. Among these people, seismic exploration 
practitioners are among the most concerned ones : in fact, their data have a poor or even 
very bad signal-to-noise ratio. This is the reason why correlated noise filtering 
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techniques are an important part of seismic data processing softwares. The most 
conventional techniques use a transform (Fourier transform for example) where signal 
and noise are located in different areas of space, thus allowing separation of the signal 
and of the noise. For a general presentation of the conventional methods for noise 
5 elimination on seismic data, one may refer to the book by Yilmaz : 

- Yilmaz, O. 1987 : « Seismic data processing », Investigation in Geophysics No.2, 
Society of Exploration Geophysicists, Tulsa, 1987. 

However, these techniques are not perfect : they presuppose that a transformation 
allowing complete separation of the signal and of the noise has been found. In this 

10 context, correlated noises are particularly bothersome because they can be difficult to 
separate from the signal (which is also correlated) and can have high amplitudes. It is 
therefore often difficult to manage the following compromise : either the signal is 
preserved, but a large noise residue remains, or the noise is eliminated, but then the 
signal is distorted. These filtering techniques can be implemented before solving the 

15 inverse problem : they constitute then a preprocessing applied to the data. The quality of 
the inverse problem solution then widely depends on the ability of the filters to 
eliminate the noise without distorting the signal. 

The approach introduced by Nemeth et al. in the following publication : 

Nemeth T. et al. (1999), « Least-Square Migration of Incomplete Reflection Data », 
20 Geophysics, 64, 208-221 

constitutes an important advance for the inversion of data disturbed by high-amplitude 
correlated noises. These authors propose eliminating the noise by solving an inverse 
problem : after defining the space of the correlated noises as the image space of a vector 
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space B (space of the noise-generating functions) by a linear operator T and the signal 
as the image space of a vector space M (models space) by a modelling operator F 
(assumed to be linear), they seek the signal in F(M) and the noise in T(B) whose sum is 
as close as possible (in the sense of the Euclidean norm or, on a continuous basis, of the 

5 norm L 2 ) to the measured datum. This technique constitutes an important advance 
insofar as it allows elimination (or at least very substantial reduction) of high-amplitude 
correlated noises (surface waves for example) that are difficult to separate from the 
signal by means of conventional filtering techniques. However, according to these 
authors, an increase by an order of magnitude of the computing time required for 

10 solution of the conventional inverse problem, i.e. without seeking the correlated 
component of the noise, is the price to pay for these performances. Furthermore, the 
result obtained with the method is extremeley sensitive to any inaccuracy introduced 
upon definition of operator T : this is the ineluctable compensation for the high aptitude 
of the method to discriminate signal and noise. 

15 SUMMARY OF THE INVENTION 

The method according to the invention allows to estimate, from data obtained by 
exploration of an underground zone, a model representative of the distribution, in the 
zone, of at least one physical quantity, this model being amply free of the presence of 
correlated noises that may be contained in the data. The method essentially comprises 
20 the following stages : 

a) acquisition of measurements giving information about certain physical 
characteristics of the zone by following a predetermined experimental protocol, 
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b) specification of a modelling operator which associates, with a model of each 
physical quantity, synthetic data that constitute the response of the model, the 
measurements and the synthetic data belonging to a data space, 

c) for each correlated noise referenced by a subscript j ranging from 1 to J, 
5 selection of a modelling operator which associates a correlated noise with a noise- 
generating function belonging to a predetermined space of noise-generating functions, 

d) specification of a norm or of a semi-norm in the data space, 

e) specification of a semi-norm in the noise-generating functions space for which 
each noise modelling operator establishes substantially an isometric relation between 

10 the noise-generating functions space and the data space, 

f) definition of a cost function quantifying the difference between the 
. measurements on the one hand and the superposition of the model response and of the 

correlated noises associated with the noise-generating functions on the other hand, and 

g) adjustment of the model and of the noise-generating functions by minimizing the 
15 cost function, by means of an algorithmic method taking advantage of the isometry 
properties of the noise modelling operators. 

According to a first implementation mode, the distribution as a function of the depth 
of the acoustic impedance in the medium is sought, the correlated noises affecting the 
data are tube waves identified each by parameters characterizing their propagation, the 
20 measured data are VSP type data obtained by means of pickups suited to detect the 
displacement of particles in the medium in response to a localized seismic excitation, 
the location of the pickups, the recording time and the time sampling points being 
defined, and the modelling operator selected associates the synthetic data with an 
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acoustic impedance distribution as a function of the evaluated depth in traveltime and 
with the vertical stress measured as a function of time at the depth of the first pickup. 

The cost function selected to quantify the difference between the measurements on 
the one hand and the superposition of the model response and of the correlated noises 
5 associated with the noise-generating functions on the other hand is for example the 
square of the semi-norm of this difference in the data space. 

According to an embodiment, adjustment of the model and of the noise-generating 
functions is for example obtained by means of a relaxation method for eliminating the 
unknowns corresponding to each correlated noise generating function, this relaxation 
10 method being implemented within the iterations of a quasi-Newtonian algorithm for 
calculation of the model. 

According to an embodiment, numerical calculation of the image of a model by this 
operator is carried out for example by numerical solution of the ID waves equation for 
the model considered, by selecting values taken by the displacement of the particles at 
15 the locations of pickups and at the previously defined time sampling points, and by 
applying an operator likely to compensate for the spherical divergence and attenuation 
effects. 

According to an embodiment, the numerical noise modelling operator is a finite- 
difference centered numerical scheme for discretizing the noise transport equation, and 
20 the noise-generating function involved as the initial condition along the edge of said 
zone associates with each correlated noise a noise-generating functions space consisting 
of the support time functions in a given time interval 
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Various examples of norms or semi-norms selected for the data space and the noise- 
generating functions space are given in the detailed description hereunder. 

According to a second implementation mode, the distribution of disturbances, in 
relation to a previously selected reference medium, of the impedance and of the velocity 

5 in said zone of the medium is sought, the correlated noises affecting the data are due to 
multiple reflections whose kinematics and amplitude variations with the offset have 
been previously estimated, the measured data are picked up by seismic surface pickups, 
the location of said pickups, the seismic excitation mode, the recording time and the 
time sampling points being defined, and the modelling operator being defined via a 

10 linearization of the waves equation around the reference medium. 

The cost function selected to quantify the difference between the measurements on 
the one hand and the superposition of the model response and of the correlated noises 
associated with the noise-generating functions on the other hand is for example the 
square of the semi-norm of this difference in the data space. 

15 According to an embodiment, adjustment of the model and of the noise-generating 

functions is obtained by means of a block relaxation method for eliminating the 
unknowns corresponding to each correlated noise generating function, this relaxation 
method being implemented within the iterations of a conjugate gradient algorithm for 
calculation of the model. 
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BRIEF DESCRIPTION OF THE FIGURES 

Other features and advantages of the method according to the invention will be 
clear from reading the description hereafter of an embodiment given by way of non 
limitative example, with reference to the accompanying drawings wherein : 

- Figure 1 shows an example of data obtained by means of a vertical seismic 
prospecting method (VSP), 

- Figure 2 shows a model of distribution of the acoustic impedance with depth, 

- Figure 3 shows synthetic VSP data obtained on the basis of the impedance model of 
Fig.2, 

- Figure 4 shows an example of VSP data contaminated by a single correlated noise, 

- Figure 5 shows an example of VSP data contaminated by two correlated noises, 

- Figure 6 shows the impedance distribution as a function of depth, obtained by 
inversion of the noise-containing data of Fig.4, 

- Figure 7 shows the inversion residues (differences between the data of Fig.4 and the 
seismic response of the model of Fig.6), 

- Figure 8 shows the impedance distribution as a function of depth, obtained by 
inversion of the noise-containing data of Fig.5, 

- Figure 9 shows the corresponding inversion residues (differences between the data 
of Fig.5 and the seismic response of the model of Fig.8), 

- Figure 10 shows the impedance distribution as a function of depth, obtained by 
inversion of the noise-containing data of Fig.4 by seeking the correlated noise in form 


of the superposition of two correlated noises having inaccurate propagation properties, 
Le. different from those of the noise appearing in the data of Fig.4, 

- Figure 1 1 shows the corresponding inversion residues (differences between the data 
of Fig.4 and the seismic response of the model of Fig. 10), 

- Figures 12 A, 12B respectively show seismic data with multiple reflections and the 
seismic response of the model obtained after conventional linearized inversion, 

Figures 13 A, 13B respectively show an example of impedance distribution and of 
propagation velocity for which the seismic data of Fig.l2A constitute the seismic 
response, and 

- Figures 14 A, 14B respectively show the comparison of the seismic responses of the 
models obtained by conventional inversion (14A) and by inversion according to the 
proposed method involving picking of the moveout of the multiple and an estimation of 
the amplitude variations with the offset (14B). 

DETAILED DESCRIPTION 

Formulation of the problem 

We have data resulting from a sampled measurement of a function. This function 
depends on several variables (space or space-time variables for example). 

These data, called d, correspond to measurements carried out to obtain information 
about a model m. They contain various noises : 

- correlated noises (or a superposition of correlated noises), 

- non correlated noises. 
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The problem is to determine quantitatively model m (or functions of this model) 
from datum d. 

We therefore select a modelling operator F (linear or not) which associates with 
model m the response of the model This operator actually models the real physical 
5 phenomenon only imperfectly. This is essentially (but not exclusively) the reason why 
correlated noises appear in the data : these noises correspond to a signal related to the 
model but the relation between them appears to be too complex to be included in 
modelling operator F, which we want to keep relatively simple for various reasons. 

Under such conditions, we model the noise by introducing one or more noise- 
10 generating spaces Bj (j ranging from 1 to J) and the associated noise modelling operators 
called Tj. In order to find model m from data d, we propose seeking this model as the 
solution to the optimization problem. 

min j c(m,fl,A,^ (i) 

where || \ D is a norm (or possibly a semi-norm) in the data space. 

15 The above formalism differs from the approach proposed by Nemeth et al. (1999) in 

that : 

- modelling operator F is not necessarily linear, 

- we are justified in considering a superposition of correlated noises of different types 
(in the sense that they are modelled by different operators referenced by subscripts j), at 

20 the price of a marked increase in the number of calculations concerning the number of 
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unknowns (which particularly complicates the solution of the optimization problem) as 
well as the number of noise modellings to be carried out, 

- we reserve the possibility of selecting norm || || 0 or possibly the semi-norm to suit 
our convenience in the data space. 

The choice that we are going to make for || \\ D will be imposed by considerations 

linked with the efficiency of the solution method (described below), this efficiency 
allowing to deal with the case of a correlated noise superposition. This is an important 
opening : besides the possibility of processing data containing correlated noises with 
complex characteristics, we can in this way overcome the difficulties encountered in 
Nemeth et al's approach, linked with the high sensitivity of the result to an inaccuracy 
introduced upon definition of the noise modelling operator : to overcome this difficulty, 
we just have to seek the noise by superposition of noises associated with approximate 
modelling operators Tj. 

It can be noted that selection of the square of norm j| |j D involved in the 

optimization problem is not essential : we do not change the solution by selecting 
instead a function which, like the square function, is an increasing function on 9T (or a 
decreasing function provided that min is changed into max). 

The solution method 

The method consists in suitably selecting : 

- norm (or semi-norm) || || D in the data space, 

- norms || J in each space Bj 
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so as to obtain a high-performance algorithmic solution to optimization problem (1). 
The expression « suitably select » means see to it that each operator Tj constitutes an 
isometry (or an isometry approximation) for norms || || and || || D respectively in the 

initial space and in the end space. The main idea in the determination of norms allowing 
to make each operator Tj isometric is based on the existence of energy conservation type 
equalities verified by the solutions to the partial differential equations (or of discrete 
energy equalities for the solutions to certain numerical schemes used for discretization 
of these equations). 

If we can make such choices, it is possible to make the cost function minimization 
very simple by eliminating the unknowns associated with the correlated noise 
generating functions (determination of Schur's complement). This elimination is 
peformed by means of a suitable algorithm such as the block relaxation method, a block 
being associated with a noise-generating function, or the conjugate gradient method 
with symmetrical-block Gauss-Seidel preconditioning, etc., all of them algorithmic 
implementation methods well-known to the man skilled in the art, presented for 
example by : 

- Golub, G.H. et aL (1983) « Resolution numerique des grands systemes lineaires» 
(Eyrolles). 

Simplified minimization of the cost function allows to considerably reduce the time 
required for numerical solution. 
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First implementiation example applied to ID inversion of VSP data con taminated hy 
a tube wave 

VSP (Vertical Seismic Profile) data belong to conventional well survey 
acquisitions. They have multiple applications, one of the main ones being their use for 
establishing a connection between seismic surface data and logging measurements. The 
inversion of VSP data as described for example by : 

- Mace, D. and Lailly, P. (1984) A Solution of an Inverse problem with ID Wave 
Equation applied to the Inversion of Vertical Seismic Profile ; Proc. of the 6 th Intern. 
Conf. In « Analysis and Optimisation of Systems », Nice, 2, 309-323 

has been developed to that end. In the ID inverse problem, we assume that the 
subsurface model is laterally invariant and that a plane-wave excitation propagates 
vertically. The unknowns of the problem are the acoustic impedance distribution as a 
function of depth (measured in vertical traveltime and over an interval starting at the 
depth of the first receiver) and the seismic excitation mode (precisely the time function 
characterizing the boundary condition at the depth of the first pickup). This is a non 
linear inverse problem, the VSP data depending non linearly on the impedance 
distribution. 

VSP data are often contaminated by a correlated noise of very high amplitude : the 
fluid wave. This wave is guided in the mud that invades the well. Its propagation 
velocity is low in relation to the propagation velocity in the rocks surrounding the well. 
Furthermore, the propagation velocity of this wave depends on the frequency 
(dispersive propagation). Finally, this wave reflects notably at the well bottom and it 
furthermore gives rise to another propagation mode. Typical VSP data are given in the 
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figure below : the main characteristics of the tube wave can be seen. Its amplitude and 
the extension of the contaminated zone make it difficult to use the signal contained in 
the many seismic samples contaminated by the fluid wave. 

We describe hereafter the implementation of the method for ID inversion of the 
5 VSPdata. 

Experimental protocol 
It consists in specifying : 

- the various depths at which the receivers are positioned in the well (assumed to be 
vertical) : in our experiments, the pickups cover the range of depths measured in one- 

10 way time [x^ = 0.05s ; = 0.2s], with a receiver every Ax=l ms (in one-way time) ; 
we thus obtain samples numbered from 0 to I, 

- the physical nature of the measurement performed by these pickups : in our 
experiments, the pickups measure the vertical displacement resulting from the wave 
propagation, 

15 - the recording time : in our experiments, the pickups measure the vibrational state 
over the time interval [0, T=1.5s], these data being sampled every At=l ms. We thus 
obtain samples numbered from 0 to N. 

We call di n the data sample recorded at the depth x^ + i Ax and at the time n At. We 
have of course : T = N At and x^ = x^ + 1 Ax. 

20 Acquisition of data in accordance with this experimental protocol 

Our experiments were carried out from synthetic data, i.e. obtained by numerical 
modelling, itself carried out in two stages : 
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Creation of synthetic data containing no noise : this creation was performed by 
numerical solution of the ID acoustic wave equation, the model being the acoustic 
impedance distribution shown in Fig.2 (the depth is measured in one-way time) and the 
excitation mode (Neumann's boundary condition at the depth 0) being a conventional 
5 Ricker wavelet. 

The synthetic VSP data resulting from this modelling are shown in Figure 3 : the 
vertical axis representing the observation time is graduated from 0 to L5 s and the 
horizontal axis representing the depth (in one-way traveltime) of the various pickups is 
graduated from 0.05 to 0.2 s. 

10 Addition of noise to the data 

We are going to disturb the VSP data thus calculated by adding thereto, on the one 
hand, a random noise (but filtered in the seismic frequency band) of rather high 
amplitude and, on the other hand, one or more correlated noises of very high amplitude. 
Each noise will propagate downwards with a low velocity which depends on the 

15 frequency (dispersive propagation) : these correlated noises are supposed to represent 
fluid waves. They have been modelled by numerical solution (use of the finite- 
difference centered scheme) of a transport equation with a constant propagation 
velocity. Dispersion of the waves has been obtained using large discretization intervals 
in the finite-difference scheme. In the case of the superposition of several correlated 

20 noises, the propagation velocities associated with each noise are different. Figures 4, 5 
show the VSP data contaminated by a single correlated noise (Fig.4) and by the 
superposition of two correlated noises (Fig.5). 
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Specification of the modelling operator 

Modelling is carried out by solving the wave equation : 

-|-(<t(x)% = 0 in Jx^ ,-h»[x [0,7-] 

Ot OX ox 

with Neumann's boundary condition : 

-^ m J^(W)=^)Wo,r] 

ca: 

and the initial conditions : 

Definition of the modelling operator first requires definition of the observation 
operator which formalizes the experimental protocol described in §1. The observation 
operator thus is the operator which associates with function y(x,t), solution to the above 
system, samples y(x i5 t n ) for x i =x min + i Ax and t D =n At, i=0,. . .,1 ; n=0,. ..,N. 

The modelling operator that we call F(m) thus is the (non linear) operator which 
associates with function pair m=(a(x),h(t)) samples y(Xi,t n ) for x i =x min + i Ax and t n =n At, 
i=0,...,I ; n=0,...,N. 

Specification of the modelling operator associated with each correlated noise 

We specify in this paragraph the procedure that can be used for modelling each 
correlated noise. Each correlated noise will be characterized by correlation directions 
specific thereto, which we assume to be known, at least approximately : these 
correlation directions are specified by means of a field of correlation vectors c y (x,/) of 

components (Cj*(x,t), c/(x,t)) which has to be defined at any point of the domain [x^, 
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XnaJXfOjT]. We can, in an equivalent way, specify the correlation directions by means 
of a pencil of correlation lines that will represent the field lines associated with Cj(x,t) . 

Such correlation lines can be obtained according to the technique described in patents 
EP-354,112 (US-4,972,383) filed by the applicant, from a pick of several phases 
5 characteristic of the noise, the information being extended to the whole domain [x^, 
XroJXtOjT] by means of conventional interpolation and/or extrapolation procedures. 
Specification of the correlation directions amounts to specifying the propagation 
velocity distribution of the noise Cj(x,t) that is related to the correlation vector by the 
cf (x,t) / 

relation c.(x,t)= j X^^xy ^ e 388111116 here that there are no wave amplitude 

10 variations during propagation of the wave (a different situation is presented for the 2 nd 
application example). Modelling of a correlated noise as specified above is based on the 
following observation : a wave propagating along correlation lines without amplitude 
change is solution to the transport equation : 

^b(x, t) . c[x, t) = 0 

15 In this context, we introduce : 

a space Bj of noise-generating functions which will be the space of support functions 
pj(t) in [t^t/™"] cz [0,T] and that we extend by 0 in the whole interval [0,T], 

- noise modelling operator Tj : 
TjiPjiOeBj-^bjix^eD 

20 solution to the transport equation 
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Vb J (x t t)£ J {x,f) = 0 in [x^ , XnaK ]x [0,T] 

and meeting the initial conditions : 
\bj(x,t = 0) = 0 
\b(x = x min ,0 = P J (t) 

It can be noted that the geometry of the pencil of correlation lines (or, which 
5 amounts to the same thing, the propagation velocity distribution) cannot be any 
geometry : the correlation lines cannot intersect, or the transport equation solution 
problem would not be correctly set. For the same reasons, the correlation lines cannot 
intercept twice the edges on which the initial conditions are specified. The latter 
consideration can lead (in order to model a correlated noise corresponding to an 
10 upgoing wave) to specify the initial condition no longer at x=x fJlin but at x=x max or even at 
an intermediate value, even if it means decomposing Cauchy's problem into two 
problems associated with the domains located on either side of this intermediate edge. 

We furthermore assume that the geometry of the correlation lines is such that (more 
general situations can be considered, but the following results then take a different 
15 form) : 

Cj(x i t) = Sj(x)XT(t) 

ou encore : cj(x, *) = ^(x) et (%(x,t) = ^ 

Finally, we assume that the pencil of correlation lines that intercept at x=x mm 
interval [t j min ,t j n,ax ] do not intercept at t=T interval [x ffiill ,x ma J. 


19 


In feet, the transport equation is solved by means of a numerical scheme. It is for 
example possible to use the conventional finite-difference centered scheme. We 
therefore select discretization intervals Ax/ and At/ (which are not necessarily equal to 
Ax and At, but which we assume to be submultiples of these quantities, even if more 
general situations can be considered) and we introduce a grid whose nodes are the 
points of coordinates (x^ + i'Ax/, n'At/) with i' and n' e N. The conventional finite- 
difference centered scheme explained hereunder allows, from the initial conditions, to 
calculate step by step the different values taken by function bj(x,t) at the various nodes 
of the grid (to simplify the notations, we have left out subscript j associated with the 
correlated noise considered). 

\h ri ' +l + b n ' +l - h n ' - h n '] + 

2 A? I «' r+1 *' J 

1 \h n ' +l ■+- h n ' - h n ' +l - fc n 'l - 0 

In the above formula, quantity a/ represents the evaluation of quantity a (a is a 
generic term) at the point of coordinates (x^ + TAx/, n'A^). 

Of course, the use of discretization intervals Ax/ and At/ sufficiently small in 
relation to the wavelength is required if we want the numerical scheme to give a precise 
approximation of the function solution to the transport equation. 

However, the use of larger discretization intervals (but still smaller than the 
wavelength) allows to model more complex propagation phenomena by making the 
wave propagation velocity dependent on the frequency : dispersive propagations such as 
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the propagation of tube waves can thus be modelled. Selection of suitable discretization 
intervals for modelling a particular mode of the tube wave can be made : 

By determining, after examining the data, the propagation velocity of this mode 
(which depends here not only on space and time, but also on the frequency) ; 
tomography techniques can be used to that end, such as those presented in the following 
publication : 

- Ernst, F. et aL (2000) ; Tomography of Dispersive Media ; J. Acoustic Soc. Am., 
108, 105-116. 

By knowing the dispersive propagation properties associated with the numerical 
scheme, properties that follow from an analysis of the numerical dispersion relation, as 
described for example in the following publication : 

- Alford, R.M. et al. (1974), Accuracy of Finite Difference Modeling of the Acoustic 
Wave Equation ; Geophysics, 6, 834-842. 

Specification of a norm or of a semi-norm in the data space 

A first choice consists in taking as the norm in the data space any discrete norm 
which is meant to be an approximation (via a quadrature formula) to the norm defined 
on a continuous basis by : 

1 

IM »-(/T*f 
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We can also use the semi-norm, which is a better choice as we shall see in the next 
paragraph : 



where u is any vector of the data space (i and n are the indices representing the sample 
5 numbers respectively in space and in time). 

Besides, other choices are also possible. 

Specification, in each space B j of noise-generating functions, of a norm for which 
each noise modelling operator 7} establishes an isometric relation (or the 
approximation of an isometric relation) between the noise-generating functions space 
10 and the data space (provided with the semi-norm defined in §5) 

We specify in this paragraph the procedure that can be used to define the norm in 
the noise-generating functions space associated with each correlated noise. This 
procedure being the same for each noise, we describe it generically and leave out, in the 
formulas, subscript j associated with the noise considered. 

15 A first choice consists in providing each noise-generating functions space B with 

any discrete norm meant to be an approximation (via a quadrature formula) to the norm 
defined on a continuous basis by : 
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in which case, using the feet that, by means of the hypotheses set out in §4, we have 
bj(x,T)=0 for x € [Xn^x^aj, linear operator Tj performs an isometry between Bj and D. 
However, since calculation of TjPj is carried out numerically, the isometric character of 
operator T will not be perfect but only approximate, the approximation being all the 
5 better as discretization intervals Ax, At, Ax/, At/ are small. The feet that we then have 
only approximately an isometry will lead to lower performances concerning the 
algorithmic implementation presented in paragraph 8 hereafter. 

A better choice consists in providing each noise-generating functions space Bj with 
the semi-norm : 
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In this case, operator Tj rigorously performs an isometry between Bj and D, at least 
if Ui N = OVi = 0,1- L This results from the discrete energy equality : 


n=0 * i=0 n=0 * 


/=o *=o A/ ^ J s At 

an equality valid in the case (simplified for presentation's sake) where Ax' = Ax and 
15 At' = At and if ^.m = has been defined. 

Definition of the cost function 

We define the cost function by the formula as follows : 
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j f 

C(m,j3 l9 f} 2 ,.../}j)= F{m) + Y T jPj- d \ 

Id 


^ Seeking the model and the noise-generating functions minimizing the cost function, 

this search being carried out by an algorithmic method taking advantage, explicitly or 
implicitly, of the isometry properties of the noise modelling operators (see 6) 

5 The algorithmic method uses the specificity of Schur's complement, a specificity 

associated with the isometric character of operators Tj. This can be done for example by 
realizing that minimizing C(m,p 1? p 2 ,-., Pj) amounts to minimizing 
C\m) = C(/w, p x (jn\/3 2 (m\...Pj(jny) where ifi x (m% P 2 {m\.„pj(m)) minimizes 

C(m,Pi, p 2 ,..., pj) for given m. It can be noted that the Hessian associated with this 
10 quadratic form is Schur's complement. C'(m) can be minimized by means of any 
optimization method such as the BFGS version of a quasi-Newtonian method. 

Determination of (J3 1 (m) 9 j3 2 (m) 9 ...fi j (m) is very easy thanks to the isometric 
character of operators T j? notably if choices (2) and (3) have been respectively selected 


to define || \\ D and || \\ B . If there is only one type of noise (J=l) and if the isometry is 

15 perfect, determination of fi x (jn) simply requires evaluation of the gradient of the 
quadratic form. If there is a superposition of several noises (J>1), various algorithms 
can be considered to take advantage of the isometric character of operators Tj. 

(fi\( m )>fi2( m )>"-fij( m ) can f° r example be determined by means of a block relaxation 
method (Gauss-Seidel alias), a block being associated with a group of unknowns 
20 J3j (m) , a relaxation iteration simply requiring (still in the case of a perfect isometry) a 
single evaluation of the gradient of the quadratic form associated with the block 
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considered. For problems made difficult by the proximity of the correlation directions of 
the various noise types, a preconditioned conjugate gradient method can be considered, 
the preconditioning matrix being of symmetrical block (relaxation alias) Gauss-Seidel 
type. Then again, preconditioning is very easily implemented, solution of the problem 
5 associated with a block only requiring evaluation of the gradient of the quadratic form 
associated with the block considered. 

Figures 6 to 1 1 show the results obtained with three types of experiments using the 
method. 

In the case of an inversion of noise-containing data (Fig.4) both by superposition of 
10 a correlated noise and of a random noise, Fig.6 gives the impedance distribution 
obtained, and Fig.7 shows the inversion residues. 

In the case of an inversion of noise-containing data (Fig.5) both by superposition of 
two correlated noises and of a random noise, Fig.8 gives the impedance distribution 
obtained, and Fig.9 shows the inversion residues. 

15 In the case of an inversion of noise-containing data (Fig.4) both by superposition of 

a correlated noise and of a random noise where the propagation properties of the fluid 
wave are not precisely known, it is natural to seek a correlated noise model consisting 
of the superposition of two correlated noises having inaccurate propagation velocities 
but by which the true propagation velocity of the fluid wave is framed. Fig. 10 gives the 

20 impedance distribution obtained, and Fig.l 1 shows the inversion residues. 

Second example of application of the method to determination of a model by 
linearized inversion of multi-offset data contaminated by multiple reflections 


25 


This application differs from the first one mainly in that it illustrates the 
possibilities afforded by the method for taking into account noise amplitude variations 
along the correlation directions. Another difference lies in the propagation properties of 
the correlated noises, propagation of the multiple reflections being not dispersive. 

5 Linearized inversion is one of the sophisticated methods used by geophysicists to 

obtain a quantitative model of subsurface heterogeneities, this type of information being 
precious notably for reservoir characterization. The data given here are surface seismic 
data. This method is described for example in the following publication : 

- Bourgeois, A, et aL (1989) « The Linearized Seismic Inverse Problem: an 
10 Attractive Method for a Sharp Description of Strategic Traps » : 59 th Ann. Intern. Mtg. 
Soc. Expl. Geoph. Expanded Abstract, pp.973-976. 

Furthermore, it is necessary to have a reference model, consisting (in the context of 
an acoustic model) of a velocity distribution and of an acoustic impedance distribution : 
it is the model around which a linearization will be carried out in the waves equation 

15 (Born's approximation). In order to simplify the presentation, we take the ID case but 
this hypothesis is all the less essential as the practical applications mainly concern 3D 
seismic methods. In the ID context, all of the seismic information is contained in a 
single shotpoint. Besides, the reference model is also one-dimensional : it only depends 
on the depth. The unknowns of the problem are the acoustic impedance disturbance and 

20 velocity distribution as a function of depth. On the other hand, unlike the first 
application example, the seismic excitation mode (and notably the seismic wavelet) are 
assumed to be known. It is an inverse problem made linear by means of linearization. 
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Surface seismic data are often contaminated by multiple reflections which may have 
a rather high amplitude. Their kinematics, which is characterized by moveout in the ID 
context, is generally quite different from that of primary reflexions : these two types of 
waves have generally travelled through geologic layers having different propagation 
5 velocities. This is the reason why linearized inversion is in itself an efficient technique 
for attenuating multiple reflections : the kinematics of primary reflections (the events 
modelled by Born's approximation) is entirely defined by the velocity distribution in the 
reference model : the modelled events therefore cannot adjust to multiple reflections if 
the moveout thereof differs from that of the primary reflections. However, because of 
10 the stationarities at the zero offset, there is in practice a partial adjustment and the 
linearized inversion results appear to be contaminated by the presence of the correlated 
noise consisting of the multiple reflections. The seismic data of Fig.l2A have been 
modelled by taking account of the multiple reflections (there are only 3 primary 
reflections which reach the zero offset successively at the times : 500 ; 1000 ; 1650 ms, 
15 the 2 nd arrival having a particularly high amplitude). 

In the seismic response of the model obtained after conventional linearized 
inversion (see Fig.l2B), the high-amplitude multiple arriving at approximately 1500 ms 
at the zero offset has been only weakly attenuated and the model is contaminated by the 
multiple. It is then natural to try to use the possibilities afforded by Nemeth's method in 
20 order to better discriminate multiples and primaries, notably with the improvements 
provided by our method. 

We describe hereafter the implementation of the method for ID linearized inversion 
of surface seismic data. 
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Experimental protocol 
It consists in specifying : 

- the seismic excitation mode : in our experiments, the seismic source is a source 
point time-modulated by a function (seismic wavelet) denoted by w(t), 

5 - the various locations where the receivers are positioned : in our experiments, the 
pickups are arranged at the depth 10m and cover the offset interval [x olin =0, 
Xn^lSOOm], with a receiver every Ax=100m ; we thus obtain seismic traces numbered 
from 0 to I, 

- the pickups measure the pressure field resulting from the wave propagation, 

10 - the recording time : : in our experiments, the pickups measure the vibrational state 
over the time interval [0, T= 1800ms], these data being sampled every At=l ms. We thus 
obtain samples numbered from 0 to N. 

We call di n the data sample recorded for the offset x^ + i Ax and at the time n At. 
We have of course : T = N At and x^ = x^ + I Ax. 

15 Acquisition of data in accordance with this experimental protocol 

Our experiments were carried out from synthetic data : these data were obtained by 
numerical resolution (finite-difference method) of the 2D wave equation : 

— ^ - V— VP = 5{x,z - z s )w(t) in 91 x 91 * + x[o,r] 
ac dt c 

with the boundary condition : 
20 P(z=0)=0 

and the initial conditions : 
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/>(r = 0)~(r = o) = o 

an equation wherein : 

- x, z, t respectively designate the lateral coordinate, the depth and the time, 

- a and c respectively represent the acoustic impedance and propagation velocity 
distributions (functions of the depth). 

The wavelet selected is a conventional Ricker wavelet of central frequency 30 Hz 
and the impedance a and velocity c distributions are given in Figs. 13 A, 13B. The 
seismic response of this model therefore consists of the data represented in Fig.l2A. We 
call Fnl the operator which associates the seismic data with function pair (a(z),c(z)). 

Specification of the modelling operator 

We select a reference model described by function pair (a 0 (z),c 0 (z)). In the 
experiments presented hereafter, we have selected c 0 (z)=c(z) and a 0 (z)=cte (=1). We 
choose as the modelling operator the Jacobian operator of operator F^ at point 
(a 0 (z),c 0 (z)). In fact, since we have selected c 0 (z)=c(z), the only unknown is 8a(z)=a(z)- 
a 0 (z) and only the corresponding component of the Jacobian is important. 

The modelling operator, which we call F(5m), thus is the (linear) operator which 
associates with function pair 5m=(8a(z),Sc(z)) the samples 5y(Xi,t n ) which are the 
components of vector F (5m) for Xr^x,^ + i Ax and t n =n At, i=0, . . . ,1 ; n=0, . . . ,N. 

Specification of the modelling operator associated with each correlated noise 

We specify in this paragraph the procedure that can be used for modelling each 
correlated noise. Each correlated noise will be characterized by correlation directions 
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specific thereto, which we assume to be known, at least approximately : these 
correlation directions are specified by means of a field of correlation vectors c y (x,/) of 

components (Cj*(x,t), c/(x,t)) which has to be defined at any point of the domain [x^, 
x max ]X[0,T]. We can, in an equivalent way, specify the correlation directions by means 
of a pencil of correlation lines that will represent the field lines associated with c y (x,0 . 

Such correlation lines can be obtained according to the technique described in patents 
EP-354,112 (US-4,972,383) filed by the applicant, from a pick of several phases 
characteristic of the noise, the information being extended to the whole domain [x^, 
XmJX[0,T] by means of conventional interpolation and/or extrapolation procedures, as 
described in the aforementioned patent EP-354,112 filed by the applicant. For this 
second application example wherein we want to eliminate the aforementioned multiple, 
we pick this multiple (the amplitude peak for example), which enables us to define the 
arrival time variations as a function of the offset and we define a correlation line 
continuum by simple vertical translation of the picked line (under these conditions, 
vector c(x,t) does not depend on t). Specification of the correlation directions amounts 
to specifying the propagation velocity distribution of the noise Cj(x,t) that is related to 

the correlation vector by the relation c^Xjt) = 

Unlike the 1 st application example described above, we have here noise amplitude 
variations along the correlation line. We assume that we can estimate these amplitude 
variations from a measurement on the data or by means of theoretical considerations : 
this measurement defines a function g(x). 

The noise will thus be modelled by the composition of two operators : 
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A transport operator, as in the first application example ; 

An amplitude modulation, i.e. the multiplication of the solution to the transport 
equation by function g(x). 

If we use again the scheme given for the 1 st application example (and notably if we 
5 use the same hypotheses and notations), we introduce : 

- a space Bj of noise-generating functions which will be the space of support functions 
pj(t) in [t^t^] c [0,T] and that we extend by 0 in the whole interval [0,T] ; 

noise modelling operator Tj : 

Tj.PjWeBj^bjix^SiD 

solution to the transport equation 

Vft y (x f /)f,(x,f) = 0 in [x^ f x_]x[0,r] 
10 and meeting the initial conditions : 

te,(x,f = 0) = 0 
&(* = *n^0 = /* y (0 

In fact, the transport equation is solved by means of a numerical scheme. It is for 
example possible to use the conventional finite-difference centered scheme. We 
therefore select discretization intervals Ax/ and At/ (which are not necessarily equal to 
Ax and At, but which we assume to be submultiples of these quantities, even if more 
15 general situations can be considered) and we introduce a grid whose nodes are the 
points of coordinates (x^ + i'Ax/, n'At/) with i' and n' e N. The conventional finite- 
difference centered scheme explained hereunder allows, from the initial conditions, to 
calculate step by step the different values taken by function bj(x,t) at the various nodes 
of the grid (to simplify the notations, we have left out subscript j associated with the 
20 correlated noise considered). 
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2Aa:'r n 


In the above formula, quantity Or°' represents the evaluation of quantity a (a is a 
generic term) at the point of coordinates (x^ + i'Ax/, n'At/). 

Here, discretization intervals Ax/ and At/ have to be selected sufficiently small in 
5 relation to the wavelength because we want the numerical scheme to give a precise 
approximation of the function solution to the transport equation (we want no 
dispersion). 

The last stage of the noise modelling procedure consists in selecting the samples h*' 
that belong to the nodes common to the two grids and in multiplying them by g(iAx) : 
10 we thus obtain the various components of the vector (of the data space) T/Pj). 

Specification of a norm or of a semi-norm in the data space 

A first choice would consist in taking as the norm in the data space any discrete 
norm meant to be an approximation (via a quadrature formula) to the norm defined on a 
continuous basis by : 

f dx / dt— u 2 (x,t) 

15 



32 


We can also use the semi-norm, which is a better choice as we shall see in the next 
paragraph : 


where u is any vector of the data space (i and n are the indices representing the sample 
numbers respectively in space and in time). 

Besides, other choices are also possible. 

Specification, in each noise-generating functions space, of a norm for which each 
noise modelling operator establishes an isometric relation (or the approximation to an 
isometric relation) between the noise-generating functions space and the data space 
(provided with the semi-norm defined in §5) 

We specify in this paragraph the procedure that can be used to define the norm in 
the noise-generating functions space associated with each correlated noise. 

A first choice would consist in providing each noise-generating functions space Bj 
with any discrete norm meant to be an approximation (via a quadrature formula) to the 
norm defined on a continuous basis by (here again, we have left out subscript j to 
simplify the notations) : 


1 
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xmar. 


dxg 2 (x) 



in which case, using the fact that bj(x,T)=0 for x e [x^x^, linear operator Tj performs 
an isometry between Bj and D. However, since calculation of Tjpj is carried out 
numerically, the isometric character of operator T, will not be perfect but only 
5 approximate, the approximation being all the better as discretization intervals Ax, At, 
Ax/, At/ are small. The fact that we then have only approximately an isometry will lead 
to lower performances concerning the algorithmic implementation presented in 
paragraph 8. 

A better choice consists in providing each noise-generating functions space Bj with 
10 the semi-norm (here again, we have left out subscript j to simplify the notations) : 


In this case, operator Tj rigorously performs an isometry between B, and D, at least 
if the solution to the discretization scheme of the transport equation is identically zero 
for n=N. 
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Definition of the cost function 


We define the cost function by the formula as follows : 


34 


II f-i Wd 

It can be noted that, in the above function, F^mo) is a constant vector. 

Seeking the model and the noise-generating functions minimizing the cost Junction, 
this search being carried out by an algorithmic method taking advantage, explicitly or 
implicitly, of the isometry properties of the noise modelling operators (see §6) 

The algorithmic method uses the specificity of Schur's complement, a specificity 
associated with the isometric character of operators Tj. This can be done for example by 
realizing that minimizing C(8m,pi, p 2 ,. . •> Pj) amounts to minimizing 

C\8m) = C(Sm, p x (<5m), /? 2 (^),.../? y (<5w)) 

where (p l (Sm) 9 /J 2 (Sn) 9 .^j3 J (3n)) minimizes CCSm^, p 2 ,..., p ; ) for given 8m. It can be 

noted that the Hessian associated with this quadratic form is Schur's complement. 
C'(Sm) can be minimized by means of any optimization method ; cost function C'(8m) 
being quadratic, the use of a conjugate gradient method is here particularly appropriate. 

Determination of (J3 l (7n) 9 j3 2 (m) 9 ...fi j (m) is very easy thanks to the isometric 
character of operators T j9 notably if choices (4) and (5) have been respectively selected 
to define || fl^ and || || B . If there is only one type of noise (J=l), which is the case for our 

illustration, and if the isometry is perfect (which is also the case thanks to the choices 
made), determination of fi x {rri) simply requires evaluation of the gradient of the 
quadratic form. If there is a superposition of several noises (J>1), various algorithms 
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can be considered to take advantage of the isometric character of operators Tj as 
explained for the first application example. 

Figure 14B shows the results obtained (seismic response of the solution model) by 
implementing the method when function g(x) is selected in a very simple form (affine 
function). The improvement in relation to the conventional inversion result can be 
observed (Fig. 14 A). 

Implementation examples where the physical parameter whose subsurface 
distribution is modelled is the acoustic impedance have been described. It is clear that 
the method, in its most general definition, can be applied for seeking the distribution of 
physical quantities affected by correlated noises in any heterogeneous medium. 


